unit diffract;

interface
uses common,ppot,d04_header;
Const
   nfitactions = 20;
Var
   sigmadr,lambdafit,etarg : myreal;
   chimin		   : myreal;
   etargs,lambdas	   : array[1..tloops] of myreal;
   nfita		   : integer;
   fita			   : array[1..nfitactions] of i2;
   setapsig		   : integer;

procedure dfit_init;
procedure read_dfit_inp(Var dmcs,cmcs,hloops,nte :integer);
procedure dcmc(nsteps,kk : integer; Var epot,ueff,rf,rwf,mychi:myreal);
procedure dfit(xyzfit,dwfit : integer; Var epot,rf,rwf:myreal;Var chi3:v3);
procedure randomize_params(what	:integer);
procedure dfit_start_info(var txpars : text);
procedure dfit_out(fname,msg : mystring ; epot,ueff,rwf,rf :myreal;chi3:v3);
procedure dpars_out(msg :mystring ; epot,ueff,chin,rwf,rf :myreal);
procedure paradump(var tx : text;ii,jj:integer;epot,ueff,rf,rwf:myreal;chi3:v3);
procedure hkl_out(fn:mystring);
procedure internal_rwfac(Var rfac,rwfac,mchi:myreal);
	    
implementation
{uses ppot,d04_header;}
uses tools,etools,relax,outp;

CONST
   {total number of parameters}
   map  = 6*t_ele*t_orb;
   tndw = 10; {different types of DWs}
TYPE
   glmma    = ARRAY [1..map] OF MYREAL;
   rorb_ch  = ARRAY [1..t_orb,1..t_ele] of myreal;
   r3orb_ch = ARRAY [1..t_orb,1..t_ele] of v3;
CONST
   appnum  = 40; {number of the independent approx. HKs making deca orbit}
   ndatap  = 3000;
   edata   = ndatap;
   mp	   = 1; {??}
   cromdim = z_ele;
TYPE
   v6 = array[1..6] of myreal;
   {GAUSSJ, LFIT, MRQMIN, MRQCOF, COVRST...}
   glnpbynp   = ARRAY [1..map,1..map] OF MYREAL;
   {glnpbymp  = ARRAY [1..map,1..mp] OF MYREAL;}
   glnpbymp   = glnpbynp;
   glnp	      = ARRAY [1..map] OF integer;
   gllista    = glnp;
   glcovar    = glnpbynp;
   glndata    = ARRAY [1..ndatap] OF MYREAL;
   glncabynca = glnpbymp;
   glnalbynal = glnpbynp;
{MRQMIN}
VAR
   glochisq : MYREAL;
   glbeta   : glmma;
{MY VARIABLES}
VAR
   alamu,dchi			: real;
   plist,st,st_ini,st_fix	: gllista;
   npfix			: integer;
   ppar,ppard,ppar0,alim0,alim1	: glmma;
   covar,alpha			: glcovar;
   maorb,exco,excro		: INTEGER;
   force_ft_init		: boolean;
   
   sf : ARRAY[1..cromdim,1..9] of MYREAL;
   hkl : ARRAY[1..ndatap] of i3;
   hkls: array[1..ndatap,1..appnum] of i3;
   napp,phkl: ARRAY[1..ndatap] of INTEGER;
   fo,fos,fosd,fc,fcs:glndata; {f-obs, f-obs-sigma, f-calc}
   q,qpe,dq,qp2,qd2,qm2,ex,qpe2 : glndata;
   qxyz:array[1..ndatap] of v3;
   ft:array[1..ndatap,1..t_orb,1..t_ele] of v2;
   ff: array[1..ndatap,1..t_orb,1..t_ele] of myreal;
   ftch:array[1..ndatap,1..t_ele] of v2;
   p2o:array[1..map,1..3] of integer;
   p2p:array[1..t_orb*3,1..t_ele] of integer;
   o2p:array[1..t_orb,1..t_ele,1..tndw] of integer;
   cro : ARRAY[1..ndatap,1..z_ele] of MYREAL;
   ppo:array[1..t_ele*t_orb*3,1..t_ele] of integer; {pointer to actual parameter #}
   dppo:array[1..map] of v2;
   {limits}
   psca0,psca1 : myreal; {scale factor}
   dwi0,dwi1 : myreal; {iso DW}
   dwa0,dwa1 : myreal; {aniso DW}
   dwp0,dwp1 : myreal; {perp DW}
   ac0,ac1   : myreal; {additive constant}
   adwp,awt:array[1..t_orb,1..z_ele] of myreal;
   adwa:array[1..t_orb,1..z_ele,1..tndw] of myreal;
   iadwp:array[1..t_orb,1..z_ele] of integer;
   iadwa:array[1..t_orb,1..z_ele,1..tndw] of integer;
   nelim      : integer;
   dok	      : array[1..ndatap] of boolean;
   unidw,glodw    : boolean;
   load_dw    : integer;
   lvmiters	: integer;
   rectime:myreal;
   emin,echimin	: myreal;
VAR
   adwglo		     : ARRAY[1..z_ele,1..tndw] of myreal;
   dwpg,aco,scal	     : myreal;
   qplim,sigmalim,flim	     : myreal;
   ppos,pdw,psca	     : integer;
   aniso,nndw,dwperp,aconst  : integer;
   ndwdat,nedatap	     : integer;
   ma,mfit,mfitp,mfitdw,ndat : integer;
   {this is number of orbits for each chemistry}
   dwsig:array[1..z_ele] of myreal;
   neledw:integer;
   peledw:array[1..z_ele] of integer;
   eledw:array[1..t_ele] of integer;
   pou:array[1..t_orb*t_ele*tndw] of integer;
   poe:array[1..t_orb*t_ele,1..2] of integer;
   setdwglo:boolean;
   ddatname:mystring;
   chinorm:myreal;
   nte_dfit:integer;
   apsig:myreal;
Var
   dwupdmode : integer;

function aino(i,j : integer):integer; begin aino:=naoch[i,oche[i,j]]; end;

procedure read_dfit_inp(Var dmcs,cmcs,hloops,nte :integer);
Var
   txt,txe    : text;
   o,i,j,nofr : integer;
   cmp : myreal;
begin
   emin:=1e10; echimin:=1e10;
   for i:=1 to atps do for j:=1 to tndw do adwglo[i,j]:=0;
   scal:=1.0; aco:=0.0;
   for i:=1 to z_ele do peledw[i]:=0;
   for o:=1 to t_orb do o_rel[o]:=false;
   reset(txt,'INPUT/dfit.inp');
   read(txt,nfita);
   for i:=1 to nfita do read(txt,fita[i,1],fita[i,2]); readln(txt);
   readln(txt,nte,hloops);
   nte_dfit:=nte;
   for i:=1 to nte do read(txt,lambdas[i]); readln(txt);
   for i:=1 to nte do read(txt,etargs[i]); readln(txt);
   for i:=1 to nte do read(txt,beta[i]); readln(txt);
   readln(txt,dmcs,cmcs,dwupdmode);
   readln(txt,sigmadr);
   readln(txt,cueff);
   readln(txt,o); glodw:=(o=1);
   readln(txt,aniso);
   readln(txt,dwperp);
   readln(txt,aconst);
   readln(txt,sigmalim);
   readln(txt,flim);
   readln(txt,qplim);
   readln(txt,dwi0,dwi1);
   readln(txt,dwa0,dwa1);
   readln(txt,dwp0,dwp1);
   readln(txt,ac0,ac1);
   for i:=1 to z_ele do dwsig[i]:=1.0;
   unidw:=false; setdwglo:=false;
   read(txt,i); if(i>0) then begin
      unidw:=true;
      setdwglo:=true;
      read(txt,cmp);
      for i:=1 to z_ele do dwsig[i]:=cmp;
   end; readln(txt);
   read(txt,i);{reserved for dwdatapoints sigmas}
   if(i>0) then begin
      unidw:=true;
      for o:=1 to i do begin
	 read(txt,j);
	 read(txt,dwsig[j]);
	 peledw[j]:=o;
	 eledw[o]:=j;
      end;
   end;
   readln(txt);
   neledw:=i;
   readln(txt,lvmiters,alamu,dchi); {iter alamu chi}
   read(txt,nofr);
   if (nofr>0) then for o:=1 to nofr do begin
      read(txt,i); o_rel[i]:=true;
   end;
   readln(txt);
   readln(txt,load_dw);
   readln(txt,nelim);
   readln(txt,setapsig,apsig);
   readln(txt,rcplus);
   for i:=1 to ndatap do dok[i]:=true;
   if(nelim>0) then begin
      reset(txe,'INPUT/DFIT/hkl.elim');
      for i:=1 to nelim do begin
	 readln(txe,j);
	 dok[j]:=false;
      end;
      close(txe);
   end;
   readln(txt,cromername);
   readln(txt,i); if(i=1) then debug[13]:=true; {INFO}
   close(txt);
   if (debug[13]) then begin
      writeln('aniso-dw              : ',aniso:1);
      writeln('perp-dw               : ',dwperp:1);
      writeln('aconst                : ',aconst:1);
      if (nofr=0) then begin
	 writeln('All orbits refined.');
      end else begin
	 write(' Orbits refined : ');
	 for o:=1 to nofr do if o_rel[o] then write(o:1,' ');
	 writeln;
      end;
   end;
   if (nofr=0) then for o:=1 to t_orb do o_rel[o]:=true;
   chimin:=1e10;
end;

procedure read_data(qplim,sigmalim,flim:myreal);
var
   o,i,j,k,ndt	   : integer;
   txt		   : text;
   h5		   : v5;
   fmax,pi2,qmix2  : myreal;
   qdeca2,qper2,dz : myreal;
   rb		   : array[1..2] of v2;
   w		   : v2;
   dk,qk	   : v3;
   p1,p2,p3	   : integer;
begin		
   pi2:=8.0*arctan(1.0);
   fmax:=0.0;
   reset(txt,'INPUT/DFIT/dat.hkl');
   readln(txt); {tiling name, read already}
   readln(txt,ddatname);
   {readln(txt,ndt,dk[1],dk[2],dk[3]);}
   readln(txt,ndt,rb[1,1],rb[1,2],rb[2,1],rb[2,2],dz);
   {if dblz>0 then dk[2]:=dk[2]/2.0;}
   ndat:=0;
   for i:=1 to ndt do begin
      ndat:=ndat+1;
      if (ndat>ndatap) then begin
	 writeln('Insufficient array dim ndatap=',ndatap:1);
	 halt;
      end;
      read(txt,fo[ndat],fos[ndat]);
      if (fo[ndat]>fmax) then fmax:=fo[ndat];
      {for o:=1 to 3 do begin}
	 read(txt,p1,p2,p3);
	 hkl[ndat,1]:=p1; hkl[ndat,2]:=p2; hkl[ndat,3]:=p3;
      {end;}
      for o:=1 to 3 do hkls[ndat,1,o]:=hkl[ndat,o];
      for o:=1 to 5 do read(txt,h5[o]);
      read(txt,q[ndat],qpe[ndat],dq[ndat]);
      ex[ndat]:=-q[ndat]*q[ndat]/4.0; {iso real-space}
      qpe2[ndat]:=-qpe[ndat]*qpe[ndat]/4.0; {perp-space}
      {for o:=1 to 3 do qk[o]:=hkl[ndat,o]*dk[o];}
      qk[1]:=hkl[ndat,1]*rb[1,1]+hkl[ndat,2]*rb[2,1];
      qk[2]:=hkl[ndat,1]*rb[1,2]+hkl[ndat,2]*rb[2,2];
      qk[3]:=dz;
      qxyz[ndat]:=qk;
      {aniso, nndw...}
      nndw:=1;
      qdeca2:=-(qk[1]*qk[1]+qk[2]*qk[2])/4.0;
      qper2:=-qk[3]*qk[3]/4.0;
      qmix2:=-sqrt(qk[1]*qk[1]+qk[2]*qk[2])*qk[3]/2.0;
      {aniso=1}
      if (aniso=1) then begin
	 qd2[ndat]:=qdeca2;
	 qp2[ndat]:=qper2;
	 nndw:=2;
      end;
      {aniso=2: after Yamamoto}
      if (aniso=2) then begin
	 qd2[ndat]:=qdeca2+qper2;
	 qp2[ndat]:=qdeca2-qper2;
	 nndw:=2;
      end;
      if (aniso=3) then begin
	 qd2[ndat]:=qdeca2;
	 qp2[ndat]:=qper2;
	 qm2[ndat]:=qmix2;
	 nndw:=3;
      end;
      read(txt,napp[i]);
      for j:=1 to napp[i] do begin
	 for o:=1 to 3 do read(txt,hkls[ndat,j+1,o]);
      end;
      readln(txt);
      {if (dblz>0) then hkl[i,2]:=hkl[i,2]*2;}
   end;
   {writeln(ndat:1,' ',fo[ndat]:1:4,' ',fos[ndat]:1:4,' ',qd2[ndat]:1:4,' ',qp2[ndat]:1:4);}
   close(txt);
   {if (debug[13]) then writeln(ndt:1,' expe data for ',ndat:1,' approximant hkls...');}
   for i:=1 to ndat do begin
      fo[i]:=fo[i]/fmax;
      fos[i]:=fos[i]/fmax;
   end;
   j:=0;
   for i:=1 to ndat do if ((qpe[i]<qplim) and (fo[i]/fos[i]>sigmalim) and (fo[i]>flim) and dok[i]) then begin
      j:=j+1;
      hkl[j]:=hkl[i];
      napp[j]:=napp[i];
      for k:=1 to napp[j] do hkls[j,k]:=hkls[i,k];
      fo[j]:=fo[i];
      fos[j]:=fos[i];
      qpe[j]:=qpe[i];
      qpe2[j]:=qpe2[i];
      q[j]:=q[i];
      qd2[j]:=qd2[i];
      qp2[j]:=qp2[i];
      qm2[j]:=qm2[i];
      qxyz[j]:=qxyz[i];
      dq[j]:=dq[i];
      ex[j]:=ex[i];
      phkl[j]:=i;
   end;
   ndat:=j;
   if (debug[13]) then 
      writeln(ndat:1,' datapoints (Q_perp<',qplim:1:1,'[A], |F|/sigma>',sigmalim:1:1,', |F|>',flim:1:4,') out of ',ndt:1);
   if (ndat=0) then begin
      writeln('0 datapoints: exiting.');
      halt;
   end;
   for i:=1 to ndat do begin fosd[i]:=fos[i]; fcs[i]:=0; end;{save}
end; { read_data }
procedure read_cromer;
var nel	: INTEGER;
   txt	: text;
begin

   if (debug[10]) then writeln('Reading ',cromername);
   reset(txt,cromername); { Cromer S-FACTOR table }
   nel:=0;
   repeat
      nel:=nel+1;
      readln(txt,sf[nel,1],sf[nel,2],sf[nel,3],sf[nel,4],
	     sf[nel,5],sf[nel,6],sf[nel,7],sf[nel,8],sf[nel,9]);
   until eof(txt);
   close(txt);
end;
function cromer(x : INTEGER;y:MYREAL):MYREAL;
{x: electron number; y: q-vector}
begin
   cromer:=sf[x,9]
   +sf[x,1]*exp(-sf[x,2]*y*y)
   +sf[x,3]*exp(-sf[x,4]*y*y)
   +sf[x,5]*exp(-sf[x,6]*y*y)
   +sf[x,7]*exp(-sf[x,8]*y*y);
end; { cromer }
procedure rwfac(fo,fc,sig : glndata;ndata:INTEGER; Var rfac,rwfac,mchi:myreal);
var
   i,j			   : INTEGER;
   dy,sum,dif,sumw,difw,wt : MYREAL;
BEGIN
{ 
RELIABILITY FACTORS (f=|F|):
R = Sum(|f_o-f_c|)/Sum(f_o)
wR = Sum(w*|f_o-f_c|)/Sum(w*f_o)
w =1/sig^2
}
   dif:=0; sum:=0;
   difw:=0; sumw:=0;
   mchi:=0;
   for i:=1 to ndata do
   begin
      sum:=sum+fo[i];
      dy:=abs(fo[i]-fc[i]);
      dif:=dif+dy;
      wt:=1/sig[i]/sig[i];
      sumw:=sumw+wt*fo[i];
      difw:=difw+wt*abs(fo[i]-fc[i]);
      mchi:=mchi+dy*dy*wt;
   end;
   rfac:=dif/sum;
   rwfac:=difw/sumw;
   {mchi:=mchi/ndata;}
   mchi:=mchi/chinorm/ndat;
END; { rwfac }
procedure chicomp(Var chi3:v3);
var
   i,p : INTEGER;
   dy,wt : MYREAL;
BEGIN
   for i:=1 to 3 do chi3[i]:=0.0;
   for i:=1 to ndat+ndwdat+nedatap do begin
      dy:=abs(fo[i]-fc[i]);
      wt:=1/fos[i]/fos[i];
      p:=1; if (i>ndat)then p:=2; if (i>ndat+ndwdat)then p:=3;
      chi3[p]:=chi3[p]+dy*dy*wt;
   end;
   for i:=1 to 3 do chi3[i]:=chi3[i]/chinorm/ndat;
END; { rwfac }
function poc(orb,ch : integer):integer;
Var
   equ:boolean;
   o,i:integer;
begin
   i:=0; repeat i:=i+1;
      equ:=(oche[orb,i]=ch);
   until ((equ) or (i=noche[orb]));
   if equ then poc:=i else begin
      write('POC: chemistry mismatch: ',orb:1,' ',ch:1,' -> ');
      for o:=1 to noche[orb] do write(oche[orb,o]:1,' ');
      writeln;
      halt;
   end;
end;
procedure ppar2xyz(x : glmma;Var k:integer);
var
   i,j: integer;
begin
   k:=0;
   for i:=1 to n_pxyz do begin
      for j:=1 to nn_pxyz[i] do begin
	 k:=k+1;
	 {writeln('ppar2xyz : ',k:2,' ',x[k]:10:5,dpxyz[i,j]:10:5);}
	 dpxyz[i,j]:=x[k];
      end;
   end;
end; { ppar2xyz }
procedure xyz2ppar(Var x : glmma);
var
   i,j,k: integer;
begin
   k:=0;
   for i:=1 to n_pxyz do begin
      for j:=1 to nn_pxyz[i] do begin
	 k:=k+1;
	 x[k]:=dpxyz[i,j];
      end;
   end;
end; { xyz2ppar }
procedure make_para_list(st:gllista) ;
Var
   i,j : integer;
begin
   {CONSTRUCT LISTA}
   mfit:=0;
   for i:=1 to ma do if (st[i]>0) then
   begin
      mfit:=mfit+1;
      plist[mfit]:=i;
      if (i<pdw) then mfitp:=mfit;
   end;
   j:=0;
   for i:=1 to ma do if (st[i]=0) then
   begin
      j:=j+1;
      plist[mfit+j]:=i;
   end;
   ma:=mfit+j;
   {writeln('make_para_list : ma ',ma:1,' mfit ',mfit:1,' ',mfitp:1);}
end; { make_para_list }
procedure para2avedw;
Var
   i,j,k,l,m,ch: integer;
begin
   for i:=1 to z_ele do for j:=1 to nndw do adwglo[i,j]:=0.0;
   k:=pdw-1;
   for i:=1 to maxorb do begin
      for j:=1 to noche[i] do begin
	 ch:=oche[i,j];
	 for l:=1 to nndw do begin
	    k:=k+1;
	    adwglo[ch,l]:=adwglo[ch,l]+ppar[k]*aino(i,j);
	 end;
	 if (dwperp=2) then begin
	    k:=k+1;
	 end;
      end;
   end;
   for i:=1 to nnele do if (nnat[ele[i]]>0) then begin
      for l:=1 to nndw do begin
	 adwglo[ele[i],l]:=adwglo[ele[i],l]/nnat[ele[i]];
      end;
   end;
end; { para2avedw }
procedure unidw_setup(Var nuni :integer );
Var
   i,j,k,l,m : integer;
begin 
   if setdwglo then begin
      neledw:=nnele;
      for i:=1 to nnele do eledw[i]:=ele[i];
      for i:=1 to z_ele do peledw[i]:=pele[i];
   end;
   for i:=1 to t_orb*t_ele*tndw do
      for j:=1 to 3 do pou[i]:=0;
   k:=0; {dw-datapoint counter}
   if(pdw>0) then begin
      m:=0; {parameter counter}
      for i:=1 to maxorb do begin
	 for j:=1 to noche[i] do begin
	    for l:=1 to nndw do begin
	       m:=m+1;
	       if((peledw[oche[i,j]]>0) and (aino(i,j)>0))then begin
		  k:=k+1;
		  pou[k]:=m;
		  fo[ndat+k]:=adwglo[oche[i,j],l];
		  fos[ndat+k]:=dwsig[oche[i,j]];
	       end;
	    end;
	 end;
      end;
   end;
   nuni:=k;
end; { unidw_setup }
procedure efit_setup(nuni :integer; Var nene :integer);
Var
   i,j,k,o : integer;
begin
   k:=0;
   for i:=1 to n_pxyz do
      for j:=1 to nn_pxyz[i] do begin
	 k:=k+1;
	 fo[ndat+nuni+k]:=0;
	 fos[ndat+nuni+k]:=lambdafit/sqrt(ndat*chinorm);
      end;
   nene:=k;
end;
procedure para2orb; {only need to transfer DW parameters}
Var
   i,j,k,l : integer;
begin
   k:=pdw-1;
   {writeln('para2orb: ',pdw:1);}
   for i:=1 to maxorb do begin
      for j:=1 to noche[i] do begin
	 for l:=1 to nndw do begin
	    k:=k+1;
	    adwa[i,oche[i,j],l]:=ppar[k];
	 end;
	 if (dwperp=2) then begin
	    k:=k+1;
	    adwp[i,oche[i,j]]:=ppar[k];
	 end;
      end;
   end;
   if (dwperp=1) then begin
      k:=k+1;
      dwpg:=ppar[k];
   end;
   if (k<>psca-1) then begin
      writeln('para2orb: mismatch in number of DW parameters.');
      halt;
   end;
   scal:=ppar[k+1];
   if (aconst=1) then aco:=ppar[k+2];
end;
procedure para2orbnew; {only need to transfer DW parameters}
Var
   i,j,k,l : integer;
begin
   if(pdw>0) then begin
      k:=pdw-1;
      {writeln('para2orb: ',pdw:1);}
      for i:=1 to maxorb do begin
	 for j:=1 to noche[i] do begin
	    for l:=1 to nndw do begin
	       k:=k+1;
	       adwa[i,oche[i,j],l]:=ppar[k];
	    end;
	    if (dwperp=2) then begin
	       k:=k+1;
	       adwp[i,oche[i,j]]:=ppar[k];
	    end;
	 end;
      end;
   end;
   if (dwperp=1) then dwpg:=ppar[psca-1];
   scal:=ppar[psca];
   if (aconst=1) then aco:=ppar[psca+1];
end;
procedure orb2para(move	:integer);
var
   o,i,j,k,l,myst : integer;
begin
   for i:=1 to map do
      for j:=1 to 3 do p2o[i,j]:=0;
   for i:=1 to maxorb do
      for j:=1 to noche[i] do
	 for o:=1 to 2 do o2p[i,j,o]:=0;
   ma:=0; k:=0;
   {positional parameters}
   if (move>0) then begin
      ppos:=1; mfitp:=0;
      for i:=1 to n_pxyz do begin
	 for j:=1 to nn_pxyz[i] do ppo[i,j]:=0;
	 for j:=1 to nn_pxyz[i] do  begin
	    ma:=ma+1;
	    ppar[ma]:=dpxyz[i,j];
	    ppo[i,j]:=ma;
            p2o[ma,1]:=pxyz2orb[i];
	    p2o[ma,2]:=ch_pxyz[i,j];
	    p2p[ma,1]:=i;
	    p2p[ma,2]:=j;
	    if (aino(pxyz2orb[i],j)>0) then st[ma]:=1 else st[ma]:=0;
	    if (st[ma]=1) then begin k:=k+1; end;
	    alim0[ma]:=-1.0;
	    alim1[ma]:=1.0; 
	 end;
      end;
   end else ppos:=0;
   {DW parameters}
   pdw:=ma+1; k:=0; mfitdw:=0;
   for i:=1 to maxorb do begin
      for j:=1 to noche[i] do begin
	 myst:=aino(i,j);
	 for l:=1 to nndw do begin
	    ma:=ma+1;
	    p2o[ma,1]:=i; p2o[ma,2]:=oche[i,j]; p2o[ma,3]:=l;
	    o2p[i,j,l]:=ma;
	    ppar[ma]:=adwa[i,oche[i,j],l];
	    st[ma]:=iadwa[i,oche[i,j],l];
	    if (myst=0) then st[ma]:=0;
	    if (st[ma]=1) then k:=k+1;
	    alim0[ma]:=dwa0;
	    alim1[ma]:=dwa1;
	 end;
	 if (dwperp=2) then begin
	    ma:=ma+1;
	    p2o[ma,1]:=i; p2o[ma,2]:=oche[i,j]; p2o[ma,3]:=nndw+1;
	    ppar[ma]:=adwp[i,oche[i,j]];
	    st[ma]:=iadwp[i,oche[i,j]];
	    if (myst=0) then st[ma]:=0;
	    if (st[ma]=1) then begin k:=k+1; mfitdw:=mfitdw+1; end;
	    alim0[ma]:=dwp0;
	    alim1[ma]:=dwp1;
	 end;
      end;
   end;
   if (dwperp=1) then begin
      ma:=ma+1;
      p2o[ma,1]:=0; p2o[ma,2]:=2;
      ppar[ma]:=dwpg;
      st[ma]:=1;
      alim0[ma]:=dwp0;
      alim1[ma]:=dwp1;
   end;
   {scale factor}
   ma:=ma+1;
   p2o[ma,1]:=0; p2o[ma,2]:=1;
   psca:=ma;
   ppar[ma]:=scal;
   st[ma]:=1;
   alim0[ma]:=0;
   alim1[ma]:=10000;
   {constant F_calc = sqrt(F_obs^2+aconst)}
   if (aconst=1) then begin
      ma:=ma+1;
      p2o[ma,1]:=0; p2o[ma,2]:=3;
      ppar[ma]:=aco;
      st[ma]:=1;
      alim0[ma]:=ac0;
      alim1[ma]:=ac1;
   end;
end; { orb2para }
procedure glo2para;
var
   i : integer;
begin
   ma:=0;
   {DW parameters}
   for i:=1 to z_ele do if (nnat[i]>0) then begin
      ma:=ma+1;
      ppar[ma]:=adwglo[i,1];
      alim0[ma]:=dwa0;
      alim1[ma]:=dwa1;
      st[ma]:=1;
   end;
   if(aniso>0) then begin
      for i:=1 to z_ele do if (nnat[i]>0) then begin
	 ma:=ma+1;
	 ppar[ma]:=adwglo[i,2];
	 alim0[ma]:=dwa0;
	 alim1[ma]:=dwa1;
	 st[ma]:=1;
      end;
   end;
   if (dwperp=1) then begin
      ma:=ma+1;
      ppar[ma]:=0.0;
      st[ma]:=1;
      alim0[ma]:=dwp0;
      alim1[ma]:=dwp1;
   end;
   {scale factor}
   ma:=ma+1;
   psca:=ma;
   ppar[ma]:=scal;
   st[ma]:=1;
   alim0[ma]:=0;
   alim1[ma]:=10000;
   for i:=1 to ma do plist[i]:=i;
   mfit:=ma;
end; { glo2para }
procedure para2glo;
var
   i,k : integer;
begin
   k:=0;
   {global DW parameters}
   for i:=1 to z_ele do if (nnat[i]>0) then begin
      k:=k+1;
      adwglo[i,1]:=ppar[k];
      {writeln('ISO para2glo: ',k:1,' ',ppar[k]:1:5);}
   end;
   if (nndw>1) then begin
      for i:=1 to z_ele do if (nnat[i]>0) then begin
	 k:=k+1;
	 adwglo[i,2]:=ppar[k];
	 {writeln('ANISO para2glo: ',k:1,' ',ppar[k]:1:5);}
      end;
   end;
   if (dwperp=1) then begin
      k:=k+1;
      dwpg:=ppar[k];
      {writeln('DWPERP para2glo: ',k:1,' ',ppar[k]:1:5);}
   end;
   {scale factor}
   k:=k+1;
   scal:=ppar[k];
end; { para2glo }
procedure funcs_mrq_elas(po: integer; a:glmma; Var yfit:myreal; Var dyda:glmma; mma:integer);
{energy datapoint: sets positional derivatives from the pair potentials}
Var
   p1,p2,i,pp,k0 : integer;
   se,sel	 : rsites;
   els,ep	 : myreal;
BEGIN
   pp:=po-ndat-ndwdat;
   yfit:=0;
   for i:=1 to mma do dyda[i]:=0;
   if(pp=1) then begin
      ppar2xyz(a,k0);
      elastic_E(cueff,ep,els,se,sel); {sets d_pxyz}
   end;
   p1:=p2p[pp,1];
   p2:=p2p[pp,2];
   if (ndps[p1,p2]>0) then begin
      yfit:=d_pxyz[p1,p2]/ndps[p1,p2];
      dyda[pp]:=-1.0/ndps[p1,p2];
   end;
END; { funcs_mrq_ene_orb }
procedure funcs_mrq_ene_one(po: integer; a:glmma; Var yfit:myreal; Var dyda:glmma; mma:integer);
{energy datapoint: sets positional derivatives from the pair potentials}
Var
   i,j,k,k0,l,o	: integer;
   se		: rsites;
   dse		: rsites3;
   na,ch,op,pch	: integer;
   epot		: myreal;
   br		: v33;
   dsp		: v3;
   dd		: glmma;
BEGIN
   for i:=1 to mma do dyda[i]:=0.0;
   if ((mma>map) or (mma<0)) then begin
      writeln('f-mrq-e : mma>map : ',mma:1,'>',map:1); halt;
   end;
   for i:=1 to mma do dd[i]:=0.0;
   PPAR2XYZ(a,k0);
   UPPW(cueff,epot,se,dse);
   na:=0;
   for i:=1 to tnodes do begin
      for j:=1 to nan[i] do if (n_a[i,j,1]>0) then begin
	 na:=na+1;
	 get_atom_basis(n_para[i,j],br);
	 mxle3(br,dse[na],dsp); {cart to parameter basis}
	 ch:=n_a[i,j,1];
	 op:=n_a[i,j,2];
	 pch:=poc(op,ch);
	 for o:=1 to 3 do if (n_para[i,j,o]>0) then begin
	    l:=ppo[n_para[i,j,o],pch]; {this is parameter number}
	    if(l>k0)then begin
	       writeln(i:1,' ',j:1,' ','funcs_mrq_ene : param pointer l=',l:1,' > k0=',k0:1);
	    end;
	    if (l>0) then dd[l]:=dd[l]+dsp[o];
	 end;
      end;
   end;
   for i:=1 to mma do begin
      if (i<=k0) then begin
	 {writeln('m_e-one : ',i:1,' ',pxyz2orb[i]:1,' ',a[i]:1:4,' ',dd[i]/na:1:4);}
	 dyda[i]:=-dd[i]/na; {why minus here?}
	 {dyda[i]:=-dd[i]/aino(p2o[i,1],p2o[i,2])/na;}
      end else dyda[i]:=0.0;
   end;
   yfit:=epot/na;
   {writeln('f-mrq-e-one : ',po:1,' ',fo[po]:1:5,' ',fos[po]:1:7,' ',yfit/na:1:7);}
   {writeln('mrq_ene: ',mma:1,' ',po:1,' ',yfit:1:4,' ',na:1);}
END; { funcs_mrq_ene }
procedure funcs_mrq_dd(po: integer; a:glmma; Var yfit:myreal; Var dyda:glmma; mma:integer);
{
po: hkl#
a: parameters
dyda: derivatives
}
var
   i,j,o,io,jo,k,l,pch : INTEGER;
   ip,ipd,ch,op,na     : integer;
   s		       : i3;
   co,si,arg,pi2,ac    : myreal;
   c,y,dwexpg,cmp      : myreal;
   x,z1,z2	       : v2;
   d		       : array[1..2] of glmma;
   dwex,dwexp,cro      : array[1..t_orb,1..t_ele] of myreal;
   reinit	       : boolean;
BEGIN
   if (mma>map) then begin
      writeln('f_m_dd : mma>map : ',mma:1,'>',map:1); halt;
   end;
   {writeln('fm-dd ',po:1);}
   for i:=1 to mma do dyda[i]:=0.0;
   pi2:=8.0*arctan(1.0);
   {set my params}
   k:=0;
   {optionally reinit structure factor (and positional derivatives)}
   reinit:=(ppos=1);
   if (reinit) then PPAR2XYZ(a,k);
   {for i:=1 to pdw-1 do if (st[i]=1) then reinit:=true;}
   {if (po=1) then writeln('funcs_mrq ',ppos:1,' ',reinit);}
   if (reinit) or (force_ft_init) then begin
      for i:=1 to maxorb do
	 for j:=1 to t_ele do
	    for o:=1 to 2 do ft[po,i,j,o]:=0.0;
      for j:=1 to pdw-1 do
	 for o:=1 to 2 do dppo[j,o]:=0.0;
      na:=0;
      for i:=1 to tnodes do begin
	 for j:=1 to nan[i] do if (n_a[i,j,1]>0) then begin
	    na:=na+1;
	    ch:=n_a[i,j,1];
	    op:=n_a[i,j,2];
	    pch:=poc(op,ch);
	    ip:=1; ipd:=1;
	    para2xyz(n_para[i,j],ch,ip,ipd,n_r[i,j]); {n_r is cryst basis}
	    s[1]:=hkl[po,1]; s[2]:=hkl[po,2]; s[3]:=hkl[po,3];
	    arg:=0.0; for o:=1 to 3 do arg:=arg+s[o]*(rt[i,o]+n_r[i,j,o]);
	    co:=cos(pi2*arg);
	    si:=sin(pi2*arg);
	    ft[po,op,pch,1]:=ft[po,op,pch,1]+co;
	    ft[po,op,pch,2]:=ft[po,op,pch,2]+si;
	    if (reinit and (mma>0)) then begin {derivatives}
	       for o:=1 to 3 do if (n_para[i,j,o]>0) then begin
		  jo:=n_para[i,j,o];
		  l:=ppo[jo,pch];
		  {if (po=5) then writeln(na:1,' ',ch:1,' ',i:1,' ',j:1,' ',jo:1,' ',l:1);}
		  if (l>0) then begin
		     if (o<3) then begin
			for io:=1 to 2 do z1[io]:=re[n_para[i,j,3+o],io];
			tra(z1,z2); {derivative in cryst basis}
			for io:=1 to 2 do
			   dppo[l,1]:=dppo[l,1]-pi2*z2[io]*s[io]*si;
			for io:=1 to 2 do
			   dppo[l,2]:=dppo[l,2]+pi2*z2[io]*s[io]*co;
		     end else begin
			dppo[l,1]:=dppo[l,1]-pi2*n_para[i,j,6]*s[3]*co*multi;
			dppo[l,2]:=dppo[l,2]+pi2*n_para[i,j,6]*s[3]*si*multi;
		     end;
		  end;
	       end;
	    end;
	 end;
      end;
      {normalization per atom}
      for i:=1 to maxorb do
	 for j:=1 to t_ele do
	    for o:=1 to 2 do ft[po,i,j,o]:=ft[po,i,j,o]/na;
      if(mma>0) then begin
	 for i:=1 to pdw-1 do begin
	    for o:=1 to 2 do
	       dppo[i,o]:=dppo[i,o]/na;
	    {if (i=24) then writeln(po:1,' ',i:1,' ',dppo[i,1]:1:5,' ',dppo[i,2]:1:5,'   ',s[1]:1,' ',s[2]:1,' ',s[3]:1);}
	 end;
      end;
   end;
   {intensity and partial derivatives}
   for o:=1 to 2 do x[o]:=0.0;
   dwexpg:=1.0;
   if (dwperp=1) then begin
      dwexpg:=exp(a[psca-1]*qpe2[po]);
   end;
   for i:=1 to maxorb do begin
      for j:=1 to noche[i] do begin
	 if (nndw=1) then begin
	    k:=k+1;
	    dwex[i,j]:=exp(a[k]*ex[po]);
	 end;
	 if (nndw=2) then begin
	    k:=k+2;
	    dwex[i,j]:=exp(a[k-1]*qd2[po]+a[k]*qp2[po]);
	 end;
	 {define qx qy qz!!}
	 if (nndw=3) then begin
	    k:=k+3;
	    dwex[i,j]:=exp(a[k-2]*qd2[po]+a[k-1]*qp2[po]+a[k]*qm2[po])
	 end;
	 dwexp[i,j]:=1.0;
	 if (dwperp=2) then begin
	    dwexp[i,j]:=exp(a[k+1]*qpe2[po]);
	 end; 
	 cro[i,j]:=cromer(oche[i,j],q[po]/2.0);
	 cmp:=dwexpg*dwex[i,j]*dwexp[i,j]*cro[i,j];
	 for o:=1 to 2 do begin
	    c:=cmp*ft[po,i,j,o];
	    x[o]:=x[o]+c;
	    if(mma>0) then begin
	       if (nndw=1) then begin
		  d[o,k]:=ex[po]*c;
	       end;
	       if (nndw=2) then begin
		  d[o,k-1]:=qd2[po]*c; {deca-plane}
		  d[o,k]:=qp2[po]*c; {periodic}
	       end;
	       if (nndw=3) then begin
		  d[o,k-2]:=qd2[po]*c; {deca-plane}
		  d[o,k-1]:=qp2[po]*c; {periodic}
		  d[o,k]:=qm2[po]*c; {mix}
	       end;
	       if (dwperp=2) then d[o,k+1]:=qpe2[po]*c;
	    end;
	 end;
	 if (dwperp=2) then k:=k+1;
      end;
   end;
   if (dwperp=1) then k:=k+1;
   if ((k+1)<>psca) then begin
      writeln('funcs_mrq: map-orb-para error: psca<>k+1. psca k+1 ',psca:1,' ',(k+1):1);
      halt;
   end;
   if (aconst=1) then ac:=a[psca+1] else ac:=0.0;
   y:=sqrt(abs(x[1]*x[1]+x[2]*x[2]+ac));
   yfit:=a[psca]*y;
   if(mma>0) then begin
      k:=0;
      {derivatives: positional parameters}
      if (ppos>0) then begin
	 for i:=1 to n_pxyz do begin
	    op:=pxyz2orb[i];
	    for j:=1 to nn_pxyz[i] do begin
	       k:=k+1;
	       dyda[k]:=a[psca]*dwex[op,j]*dwexp[op,j]*cro[op,j]*(x[1]*dppo[k,1]+x[2]*dppo[k,2])/y;
	    end;
	 end;
      end;
      {derivatives: dw factors}
      for i:=1 to maxorb do begin
	 for j:=1 to noche[i] do begin
	    for l:=1 to nndw do begin
	       k:=k+1;
	       dyda[k]:=a[psca]*(x[1]*d[1,k]+x[2]*d[2,k])/y;
	    end;
	    if (dwperp=2) then begin
	       k:=k+1;
	       dyda[k]:=a[psca]*(x[1]*d[1,k]+x[2]*d[2,k])/y;
	    end;
	 end;
      end;
      if (dwperp=1) then begin
	 k:=k+1;
	 dyda[k]:=a[psca]*qpe2[po]*(x[1]*x[1]+x[2]*x[2])/y;
      end;
      dyda[psca]:=y; {scale factor}
      if (aconst=1) then dyda[psca+1]:=a[psca]/y/2.0;
   end;
END; { funcs_mrq_dd }
procedure funcs_mrq_glo(po: integer; a:glmma; Var yfit:myreal; Var dyda:glmma; mma:integer);
{
po: hkl#
a: parameters
dyda: derivatives
}
var
   i,j,o,k,l	    : INTEGER;
   ip,ipd,ch,na	    : integer;
   s		    : i3;
   co,si,arg,pi2,ac : myreal;
   c,y,dwexpg	    : myreal;
   x		    : v2;
   d		    : array[1..2] of glmma;
   dwex,dwexp,cro   : array[1..t_ele] of myreal;
BEGIN
   for i:=1 to mma do dyda[i]:=0.0;
   if (mma>map) then begin
      writeln('f_m_glo : mma>map : ',mma:1,'>',map:1); halt;
   end;
   pi2:=8.0*arctan(1.0);
   if (force_ft_init) then begin
      na:=0;
      for j:=1 to nnele do
	 for o:=1 to 2 do ftch[po,j,o]:=0.0;
      for i:=1 to tnodes do begin
	 for j:=1 to nan[i] do if (n_a[i,j,1]>0) then begin
	    na:=na+1;
	    ch:=n_a[i,j,1];
	    ip:=1; ipd:=1;
	    para2xyz(n_para[i,j],ch,ip,ipd,n_r[i,j]); {n_r is cryst basis}
	    s[1]:=hkl[po,1]; s[2]:=hkl[po,2]; s[3]:=hkl[po,3];
	    arg:=0.0; for o:=1 to 3 do arg:=arg+s[o]*(rt[i,o]+n_r[i,j,o]);
	    co:=cos(pi2*arg);
	    si:=sin(pi2*arg);
	    ftch[po,pele[ch],1]:=ftch[po,pele[ch],1]+co;
	    ftch[po,pele[ch],2]:=ftch[po,pele[ch],2]+si;
	 end;
      end;
      {normalization per atom}
      for i:=1 to nnele do
	 for o:=1 to 2 do ftch[po,i,o]:=ftch[po,i,o]/na;
   end;
   {intensity and partial derivatives}
   for o:=1 to 2 do x[o]:=0.0;
   dwexpg:=1.0;
   if (dwperp=1) then begin
      dwexpg:=exp(a[psca-1]*qpe2[po]);
   end;
   for i:=1 to nnele do begin
      if(nndw=1) then dwex[i]:=exp(a[i]*ex[po]);
      if(nndw>1) then dwex[i]:=exp(a[i]*qd2[po]+a[i+nnele]*qp2[po]);
      cro[i]:=cromer(ele[i],q[po]/2.0);
      for o:=1 to 2 do begin
	 c:=dwexpg*dwex[i]*cro[i]*ftch[po,i,o];
	 x[o]:=x[o]+c;
	 if (nndw=1) then begin
	    d[o,i]:=ex[po]*c;
	 end;
	 if (nndw>1) then begin
	    d[o,i]:=qd2[po]*c; {deca-plane}
	    d[o,i+nnele]:=qp2[po]*c; {periodic}
	 end;
      end;
   end;
   y:=sqrt(abs(x[1]*x[1]+x[2]*x[2]));
   yfit:=a[psca]*y;
   {derivatives: dw factors}
   for i:=1 to nnele do begin
      dyda[i]:=a[psca]*(x[1]*d[1,i]+x[2]*d[2,i])/y;
      if(nndw>1)then begin
	 dyda[i+nnele]:=a[psca]*(x[1]*d[1,i+nnele]+x[2]*d[2,i+nnele])/y;
      end;
   end;
   if (dwperp=1) then begin
      dyda[psca-1]:=a[psca]*qpe2[po]*(x[1]*x[1]+x[2]*x[2])/y;
   end;
   dyda[psca]:=y; {scale factor}
END; { funcs_mrq_glo }
procedure funcs_mrq_dwave(po: integer; a:glmma; Var yfit:myreal; Var dyda:glmma; mma:integer);
Var
   i,k,pp : integer;
begin
   for i:=1 to mma do dyda[i]:=0.0;
   pp:=po-ndat;
   k:=pou[pp]+pdw-1;
   yfit:=a[k]; dyda[k]:=1.0;
   {writeln('f-mrq-dwave: ',po:1,' ',pp:1,' ',pdw:1,' ',k:1,' ',p2o[k,1]:1,' ',p2o[k,2]:1,' ',p2o[k,3]:1,' ',a[k]:1:4,' ',adwglo[p2o[k,2],p2o[k,3]]:1:4,' ',fo[ndat+pp]:1:4);}
end; { funcs_mrq_dwave }
procedure funcs_mrq(po: integer; a:glmma; Var yfit:myreal; Var dyda:glmma; mma:integer);
BEGIN
   if glodw then begin
      {GLOBAL ELEMENTAL DW FACTORS ONLY}
      funcs_mrq_glo(po,a,yfit,dyda,mma);
   end else begin
      {ENERGY DATAPOINT}
      if (po>ndat+ndwdat) then begin
	 if(nedatap=1) then funcs_mrq_ene_one(po,a,yfit,dyda,mma) else funcs_mrq_elas(po,a,yfit,dyda,mma);
      end;
      {AVEDW DATAPOINTS}
      if ((po>ndat) and (ndwdat>0) and (po<=ndat+ndwdat)) then begin
	 funcs_mrq_dwave(po,a,yfit,dyda,mma);
      end;
      {DIFFRACTION DATA}
      if (po<=ndat) then funcs_mrq_dd(po,a,yfit,dyda,mma);
   end;
END; { funcs_mrq }

PROCEDURE gaussj(VAR a: glnpbynp; n,np: INTEGER;
       VAR b: glnpbymp; m,mp: INTEGER);
{ Programs using GAUSSJ must define the types
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnpbymp = ARRAY [1..np,1..mp] OF real;
   glnp = ARRAY [1..np] OF INTEGER;
in the main routine. }
VAR
   big,dum,pivinv: MYREAL;
   i,icol,irow,j,k,l,ll: INTEGER;
   indxc,indxr,ipiv: glnp;
BEGIN
   FOR j := 1 to n DO BEGIN
      ipiv[j] := 0
   END;
   FOR i := 1 to n DO BEGIN
      big := 0.0;
      FOR j := 1 to n DO BEGIN
         IF (ipiv[j] <> 1) THEN BEGIN
            FOR k := 1 to n DO BEGIN
               IF (ipiv[k] = 0) THEN BEGIN
                  IF (abs(a[j,k]) >= big) THEN BEGIN
                     big := abs(a[j,k]);
                     irow := j;
                     icol := k
                  END
               END ELSE IF (ipiv[k] > 1) THEN BEGIN
                  {writeln('pause 1 in GAUSSJ - singular matrix'); readln}
               END
            END
         END
      END;
      ipiv[icol] := ipiv[icol]+1;
      IF (irow <> icol) THEN BEGIN
         FOR l := 1 to n DO BEGIN
            dum := a[irow,l];
            a[irow,l] := a[icol,l];
            a[icol,l] := dum
         END;
         FOR l := 1 to m DO BEGIN
            dum := b[irow,l];
            b[irow,l] := b[icol,l];
            b[icol,l] := dum
         END
      END;
      indxr[i] := irow;
      indxc[i] := icol;
      IF (a[icol,icol] = 0.0) THEN BEGIN
	 writeln('pause 2 in GAUSSJ - singular matrix : ma=',plist[icol]:1,' pdw psca ',pdw:1,' ',psca:1); {readln}
	 if (plist[icol]<(psca-1)) then 
	    writeln('ORB ',p2o[plist[icol],1]:1,' CHEM ',p2o[plist[icol],2]:1) else
	 writeln('zero=',p2o[plist[icol],1]:1,'[1=sca 2=pdw 3=aconst] :  ',p2o[plist[icol],2]:1);
	 halt;
      END;
      pivinv := 1.0/a[icol,icol];
      a[icol,icol] := 1.0;
      FOR l := 1 to n DO BEGIN
         a[icol,l] := a[icol,l]*pivinv
      END;
      FOR l := 1 to m DO BEGIN
         b[icol,l] := b[icol,l]*pivinv
      END;
      FOR ll := 1 to n DO BEGIN
         IF (ll <> icol) THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 to n DO BEGIN
               a[ll,l] := a[ll,l]-a[icol,l]*dum
            END;
            FOR l := 1 to m DO BEGIN
               b[ll,l] := b[ll,l]-b[icol,l]*dum
            END
         END
      END
   END;
   FOR l := n DOWNTO 1 DO BEGIN
      IF (indxr[l] <> indxc[l]) THEN BEGIN
         FOR k := 1 to n DO BEGIN
            dum := a[k,indxr[l]];
            a[k,indxr[l]] := a[k,indxc[l]];
            a[k,indxc[l]] := dum
         END
      END
   END
END;
PROCEDURE covsrt(VAR covar: glcovar; ncvm: INTEGER; ma: INTEGER;
       lista: gllista; mfit: INTEGER);
{ Programs using routine COVSRT must define the types
TYPE
   glcovar = ARRAY [1..ncvm,1..ncvm] OF real;
   gllista = ARRAY [1..mfit] OF INTEGER;
in the calling program. }
VAR
   j,i: INTEGER;
   swap: MYREAL;
BEGIN
   FOR j := 1 to ma-1 DO BEGIN
      FOR i := j+1 to ma DO BEGIN
         covar[i,j] := 0.0
      END
   END;
   FOR i := 1 to mfit-1 DO BEGIN
      FOR j := i+1 to mfit DO BEGIN
         IF (lista[j] > lista[i])  THEN BEGIN
            covar[lista[j],lista[i]] := covar[i,j]
         END ELSE BEGIN
            covar[lista[i],lista[j]] := covar[i,j]
         END
      END
   END;
   swap := covar[1,1];
   FOR j := 1 to ma DO BEGIN
      covar[1,j] := covar[j,j];
      covar[j,j] := 0.0
   END;
   covar[lista[1],lista[1]] := swap;
   FOR j := 2 to mfit DO BEGIN
      covar[lista[j],lista[j]] := covar[1,j]
   END;
   FOR j := 2 to ma DO BEGIN
      FOR i := 1 to j-1 DO BEGIN
         covar[i,j] := covar[j,i]
      END
   END
END;
PROCEDURE mrqcof(x,y,sig: glndata; ndata: INTEGER;
       VAR a: glmma; mma: INTEGER; lista: gllista;
       mfit: INTEGER; VAR alpha: glnalbynal;
       VAR beta: glmma; nalp: INTEGER; VAR chisq: MYREAL);
VAR
   k,j,i	    : INTEGER;
   ymod,wt,sig2i,dy : MYREAL;
   dyda		    : glmma;
   da0,da1	    : myreal;
BEGIN
   FOR j := 1 to mfit DO BEGIN
      FOR k := 1 to j DO BEGIN
         alpha[j,k] := 0.0
      END;
      beta[j] := 0.0
   END;
   chisq := 0.0;
   force_ft_init:=false; {MM}
   {writeln('mrqcof ',mma:1);}
   FOR i := 1 to ndata DO BEGIN
      funcs_mrq(i,a,ymod,dyda,mma);
      sig2i := 1.0/(sig[i]*sig[i]);
      dy := y[i]-ymod;
      FOR j := 1 to mfit DO BEGIN
         wt := dyda[lista[j]]*sig2i;
         FOR k := 1 to j DO BEGIN
            alpha[j,k] := alpha[j,k]+wt*dyda[lista[k]]
         END;
         beta[j] := beta[j]+dy*wt
      END;
      chisq := chisq+dy*dy*sig2i
   END;
   {writeln('CHISQ_TRIAL: ',chisq/ndata:1:5);}
   FOR j := 2 to mfit DO BEGIN
      FOR k := 1 to j-1 DO BEGIN
         alpha[k,j] := alpha[j,k]
      END
   END
END;
PROCEDURE mrqmin(x,y,sig: glndata; ndata: INTEGER;
         VAR a: glmma; mma: INTEGER; lista: gllista;
         mfit: INTEGER; VAR covar,alpha: glncabynca;
         nca: INTEGER; VAR chisq,alamda: MYREAL);
VAR
   k,kk,j,ihit			   : INTEGER;
   atry,da			   : glmma; oneda: glncabynca;
   plus,minus,sc,eps,da0,da1 : MYREAL;
   goto99			   : boolean;
BEGIN
   eps:=1e-6; {???}
   IF (alamda < 0.0) THEN BEGIN
      kk := mfit+1;
      FOR j := 1 TO mma DO BEGIN
	 ihit := 0;
	 FOR k := 1 TO mfit DO BEGIN
	    IF (lista[k] = j) THEN ihit := ihit+1
	 END;
	 IF (ihit = 0) THEN BEGIN
	    lista[kk] := j;  kk := kk+1
	 END ELSE IF (ihit > 1) THEN BEGIN
	       writeln('pause 1 in routine MRQMIN');
	     writeln('Improper permutation in LISTA'); readln
	 END
      END;
      IF (kk <> (mma+1)) THEN BEGIN
	 writeln('pause 2 in routine MRQMIN');
	 writeln('Improper permutation in LISTA'); readln
      END;
      alamda := 0.001;
      mrqcof(x,y,sig,ndata,a,mma,lista,mfit,alpha,glbeta,nca,chisq);
      glochisq := chisq;
      FOR j := 1 TO mma DO BEGIN
	   atry[j] := a[j]
      END
   END;
   FOR j := 1 TO mfit DO BEGIN
      FOR k := 1 TO mfit DO BEGIN
	 covar[j,k] := alpha[j,k]
      END;
      covar[j,j] := alpha[j,j]*(1.0+alamda);
      oneda[j,1] := glbeta[j]
   END;
   gaussj(covar,mfit,nca,oneda,1,1);
   FOR j := 1 TO mfit DO da[j] := oneda[j,1];
   goto99:=false;
   IF (alamda < eps) THEN BEGIN
      covsrt(covar,nca,mma,lista,mfit);
      GOTO99:=true;
   END;

   if not(goto99) then
   begin
      {TREATMENT OF LIMIT CONSTRAINTS}
      for j:=1 to mfit do
      begin
	 da0:=a[lista[j]]+da[j]-alim0[lista[j]];
	 da1:=a[lista[j]]+da[j]-alim1[lista[j]];
         if (a[lista[j]]+da[j]<alim0[lista[j]]) then
	 begin
            da[j]:=alim0[lista[j]]-a[lista[j]];
	    if(st_fix[lista[j]]=0)then npfix:=npfix+1;
	    st_fix[lista[j]]:=1;
	 end;
         if (a[lista[j]]+da[j]>alim1[lista[j]]) then
	 begin
            da[j]:=alim1[lista[j]]-a[lista[j]];
	    if(st_fix[lista[j]]=0)then npfix:=npfix+1;
	    st_fix[lista[j]]:=1;
	 end;
      end;
      
      FOR j := 1 TO mfit DO BEGIN
	 atry[lista[j]] := a[lista[j]]+da[j]
      END;

      mrqcof(x,y,sig,ndata,atry,mma,lista,mfit,covar,da,nca,chisq);
      IF (chisq < glochisq) THEN BEGIN
	 alamda := alamda/alamu;  glochisq := chisq;
	 FOR j := 1 TO mfit DO BEGIN
	    FOR k := 1 TO mfit DO BEGIN
	       alpha[j,k] := covar[j,k]
	    END;
	    glbeta[j] := da[j];
	    a[lista[j]] := atry[lista[j]]
	 END
      END
      ELSE BEGIN
	 {for j:=1 to mfit do stn[lista[j]]:=-1;}
	 alamda := alamu*alamda;  chisq := glochisq
      END;
   end;
END; { mrqmin }
procedure internal_rwfac(Var rfac,rwfac,mchi:myreal);
Var
   i,j,k,o,ndt		   : integer;
   dy,sum,dif,sumw,difw,wt : MYREAL;
   fs,fs2		   : myreal;
begin
   dif:=0; sum:=0;
   difw:=0; sumw:=0;
   mchi:=0;
   ndt:=0;
   force_ft_init:=true;
   for i:=1 to ndat do begin
      fs:=0; fs2:=0;
      for j:=napp[i] downto 1 do begin {so that in the end we have the first reinited}
	 ndt:=ndt+1;
	 for o:=1 to 3 do hkl[i,o]:=hkls[i,j,o];
	 funcs_mrq(i,ppar,fc[i],ppard,0);
	 fs:=fs+fc[i];
	 fs2:=fs2+fc[i]*fc[i];
	 sum:=sum+fo[i];
	 dy:=abs(fo[i]-fc[i]);
	 dif:=dif+dy;
	 wt:=1/fosd[i]/fosd[i];
	 sumw:=sumw+wt*fo[i];
	 difw:=difw+wt*abs(fo[i]-fc[i]);
	 mchi:=mchi+dy*dy*wt;
      end;
      if (napp[i]>1) then begin
	 fs2:=fs2/napp[i];
	 fs:=fs/napp[i];
	 fcs[i]:=sqrt(fs2-fs*fs);
      end else fcs[i]:=0;
   end;
   rfac:=dif/sum;
   rwfac:=difw/sumw;
   {mchi:=mchi/ndata;}
   mchi:=mchi/chinorm/ndt;
end;
procedure dfit_start_info(var txpars : text);
Var
   i,j	   : integer;
   tfitpos : integer;
   tfitdw  : integer;
begin
   writeln(txpars,'start_paradump');
   writeln(txpars,'info ticoname ',ticoname);
   writeln(txpars,'info deconame ',deconame);
   writeln(txpars,'info ddatname ',ddatname);
   writeln(txpars,'info ndat ',ndat:1);
   writeln(txpars,'info limits(sig qp f) ',sigmalim:1:4,' ',qplim:1:4,' ',flim:1:4);
   writeln(txpars,'info atoms ',natoms:1);
   writeln(txpars,'info tnodes ',tnodes:1);
   writeln(txpars,'info aniso aconst dwperp ',aniso:1,' ',aconst:1,' ',dwperp:1);
   writeln(txpars,'info setapsig apsig ',setapsig:1,' ',apsig:1:3);
   tfitpos:=0;
   for i:=1 to n_pxyz do
      for j:=1 to nn_pxyz[i] do if (aino(pxyz2orb[i],j)>0) then tfitpos:=tfitpos+1;
   tfitdw:=0;
   for i:=1 to maxorb do
      for j:=1 to noche[i] do if (aino(i,j)>0) then begin
	 tfitdw:=tfitdw+nndw;
	 if (dwperp=2) then tfitdw:=tfitdw+1;
      end;
   writeln(txpars,'info pospar dwpar ',mfitp:1,'/',tfitpos:1,'  ',tfitdw:1);
   writeln(txpars,'start_fit_schedule');
   if((nte_dfit>0) and (nfita>0))then begin
      if(nfita>0) then writeln('FIT SCHEDULE : ');
      for i:=1 to nfita do begin
	 if(fita[i,1]<0) then begin
	    writeln('  ',i:1,' randomize positions');
	    writeln(txpars,'  ',i:1,' randomize positions');
	 end;
	 if(fita[i,1]=1) then begin
	    writeln('  ',i:1,' xyzfit 1 no-energy-bias');
	    writeln(txpars,'  ',i:1,' xyzfit 1 no-energy-bias');
	 end;
	 if(fita[i,1]=2) then begin
	    writeln('  ',i:1,' xyzfit 2 E_targ:(lambda*100) -> ');
	    write('    ');
	    for j:=1 to nte_dfit do write(etargs[j]:1:3,':',lambdas[j]*100:1:3,' ');
	    writeln;
	    writeln(txpars,'  ',i:1,' xyzfit 2 E_targ/lambda*100 -> ');
	    write(txpars,'    ');
	    for j:=1 to nte_dfit do write(txpars,etargs[j]:1:3,'/',lambdas[j]*100:1:3,' ');
	    writeln(txpars);
	 end;
	 if(fita[i,1]=3) then begin
	    write('  ',i:1,' xyzfit 3 lambda*100 -> ');
	    for j:=1 to nte_dfit do write(lambdas[j]:1:3,' ');
	    writeln;
	    write(txpars,'  ',i:1,' xyzfit 3 lambda*100 -> ');
	    for j:=1 to nte_dfit do write(txpars,lambdas[j]:1:3,' ');
	    writeln(txpars);
	 end;
	 if(fita[i,2]<0) then writeln('  ',i:1,' randomize DW factors');
	 if(fita[i,2]<0) then writeln(txpars,'  ',i:1,' randomize DW factors');
	 if(fita[i,2]>0) then begin
	    if(fita[i,2]=2) then begin
	       writeln('  ',i:1,' optimize DW factors constrain unidw ');
	       writeln(txpars,'  ',i:1,' optimize DW factors constrain unidw ');
	    end else begin
	       writeln(i:1,' optimize DW factors');
	       writeln(txpars,i:1,' optimize DW factors');
	    end;
	 end;
      end;
      writeln('----------  fit_started  -----------------');
   end;
   writeln(txpars,'end_fit_schedule');
end;
procedure dfit_out(fname,msg: mystring ; epot,ueff,rwf,rf :myreal;chi3:v3);
var
   i,j,k,l,o,orb,ch : integer;
   npp,npdw,npa	    : integer;
   ac,se,mchi	    : myreal;
   txt		    : text;
begin
   npp:=pdw-1;
   npdw:=psca-1-npp;
   npa:=psca; if (aconst=1) then npa:=npa+1;
   rewrite(txt,fname);
   for i:=1 to 2 do begin
      for o:=1 to 5 do write(txt,bai[i,o]:5); writeln(txt);
   end;
   write(txt,tnodes:1,' ',multi:1,' ');
   writeln(txt,deconame);
   export_tiling_body(txt);
   export_params(txt);
   export_tiling_atoms(txt);
   if (dwperp=1) then dwpg:=ppar[psca-1];
   ac:=0.0; if (aconst=1) then ac:=ppar[psca+1];
   writeln(txt,'start_dfit_params -----------------------------------------');
   writeln(txt,ppar[psca]:1:5,' ',st[psca]:1,' ',ac:1:5,' ',aniso:1,' ',dwperp:1,' ',dwpg:1:4,'  # scal_fact scal_st aconst aniso dwperp dwpglo');
   k:=pdw-1;
   for i:=1 to maxorb do begin
      for j:=1 to noche[i] do begin
	 write(txt,i:2,oche[i,j]:3);
	 for l:=1 to nndw do begin
	    k:=k+1;
	    write(txt,ppar[k]:10:4,st[k]:2);
	 end;
	 if (dwperp=2) then begin
	    k:=k+1;
	    write(txt,ppar[k]:10:4,st[k]:2);
	 end;
	 writeln(txt,naoch[i,oche[i,j]]/natoms:10:4);
      end;
   end;
   writeln(txt,'end_dfit_params -------------------------------------------');
   writeln(txt,'start_dfit_rep --------------------------------------------');
   writeln(txt,'Exporter: ',msg);
   writeln(txt,ticoname,' # ticoname');
   writeln(txt,deconame,' # deconame');
   writeln(txt,ddatname,' # ddatname');
   writeln(txt,natoms:1,' ',epot/natoms:1:6,' ',ueff/natoms:1:6,' # nat epot ueff [eV/atom]');
   writeln(txt,rwf:1:6,' ',rf:1:6,' ',chi3[1]:1:3,' ',chi3[2]:1:3,' ',chi3[3]:1:3,' # wR R chi_df chi_dw chi_ene');
   writeln(txt,ndat:1,' ',npp:1,' ',npdw:1,' ',npa:1,'  # ndat pospar dwpar allpar');
   writeln(txt,aniso:1,' ',dwperp:1,' ',aconst:1,'  # aniso dwperp aconst ');
   writeln(txt,sigmalim:1:4,' ',qplim:1:4,' ',flim:1:4,'  # sigmalim qplim flim');
   writeln(txt,lambdafit:1:4,' # lambdafit');
   writeln(txt,sigmadr:1:3,' # sigmadr');
   writeln(txt,cueff:1:4,' # cueff');
   writeln(txt,nndw:1,' ',ndwdat:1,' # nndw ndwdat (unidw)');
   writeln(txt,setapsig:1,' ',apsig:1:3,' # setapsig apsig');
   for i:=1 to nnele do begin
      write(txt,ele[i]:1,' ',adwglo[ele[i],1]:1:3);
      if (nndw>0) then writeln(txt,' ',adwglo[ele[i],2]:1:3,' ',dwsig[ele[i]]:1:3,'  # aveDW(CH) aveDWani(CH) dwsig') else writeln(txt,dwsig[ele[i]]:1:3,'  # aveDW(CH) dwsig');
   end;
   writeln(txt,'end_dfit_rep --------------------------------------------');
   close(txt);
end; { dfit_out }
procedure dpars_out(msg :mystring ; epot,ueff,chin,rwf,rf :myreal);
var
   i,j,k,l,o,orb,ch : integer;
   npp,npdw,npa	    : integer;
   ac,se,mchi	    : myreal;
   txt		    : text;
begin
   {PARAMETERS OUT TO DFIT/fit}
   npp:=pdw-1;
   npdw:=psca-1-npp;
   npa:=psca; if (aconst=1) then npa:=npa+1;
   rewrite(txt,'DFIT/fit');
   writeln(txt,msg);
   writeln(txt,ticoname,' # ticoname');
   writeln(txt,deconame,' # deconame');
   writeln(txt,ddatname,' # ddatname');
   writeln(txt,natoms:1,' ',epot/natoms:1:6,' ',ueff/natoms:1:6,' # nat epot ueff [eV/atom]');
   writeln(txt,rwf:1:6,' ',rf:1:6,' ',chin:1:2,' wR R chisq');
   writeln(txt,ndat:1,' ',npp:1,' ',npdw:1,' ',npa:1,'  # ndat pospar dwpar allpar');
   writeln(txt,aniso:1,' ',dwperp:1,' ',aconst:1,'  # aniso dwperp aconst ');
   writeln(txt,sigmalim:1:4,' ',qplim:1:4,' ',flim:1:4,'  # sigmalim qplim flim');
   writeln(txt,lambdafit:1:4,' # lambdafit');
   writeln(txt,sigmadr:1:3,' # sigmadr');
   writeln(txt,cueff:1:4,' # cueff');
   writeln(txt,nndw:1,' ',ndwdat:1,' # nndw ndwdat (unidw)');
   for i:=1 to nnele do begin
      write(txt,ele[i]:1,' ',adwglo[ele[i],1]:1:3);
      if (nndw>0) then writeln(txt,' ',adwglo[ele[i],2]:1:3,' ',dwsig[ele[i]]:1:3,'  # aveDW(CH) aveDWani(CH) dwsig') else writeln(txt,dwsig[ele[i]]:1:3,'  # aveDW(CH) dwsig');
   end;
   ac:=0.0; if (aconst=1) then ac:=ppar[psca+1];
   writeln(txt,'start_dfit_params -----------------------------------------');
   writeln(txt,ppar[psca]:1:5,' ',st[psca]:1,' ',ac:1:5,' ',aniso:1,' ',dwperp:1,' ',dwpg:1:4,'  # scal_fact scal_st aconst aniso dwperp dwpglo');
   k:=pdw-1;
   for i:=1 to maxorb do begin
      for j:=1 to noche[i] do begin
	 write(txt,i:2,oche[i,j]:3);
	 for l:=1 to nndw do begin
	    k:=k+1;
	    write(txt,ppar[k]:10:4,st[k]:2);
	 end;
	 if (dwperp=2) then begin
	    k:=k+1;
	    write(txt,ppar[k]:10:4,st[k]:2);
	 end;
	 writeln(txt,naoch[i,oche[i,j]]/natoms:10:4);
      end;
   end;
   writeln(txt,'end_dfit_params -------------------------------------------');
   close(txt);
end; { dpars_out }
procedure paradump(var tx : text;ii,jj:integer;epot,ueff,rf,rwf:myreal;chi3:v3);
Var
   i,j,k,l,o : integer;
   ac	   : myreal;
begin
   ac:=0; if (aconst=1) then ac:=ppar[psca+1];
   write(tx,"_newsample_ ",ii:1,' ',jj:1,' ',epot/natoms:1:6,' ',ueff/natoms:1:6,' ',rf:1:4,' ',rwf:1:4,' ');
   for o:=1 to 3 do write(tx,chi3[o]:1:3,' ');
   writeln(tx,'# L1 L2 E Ueff rf rwf chi_df chi_dw chi_e');
   writeln(tx,'start_dfit_params -----------------------------------------');
   writeln(tx,ppar[psca]:1:5,' ',st[psca]:1,' ',ac:1:5,' ',aniso:1,' ',dwperp:1,' ',dwpg:1:4,'  # scal_fact scal_st aconst aniso dwperp dwpglo');
   k:=pdw-1;
   for i:=1 to maxorb do begin
      for j:=1 to noche[i] do begin
	 write(tx,i:2,oche[i,j]:3);
	 for l:=1 to nndw do begin
	    k:=k+1;
	    write(tx,ppar[k]:10:4,st[k]:2);
	 end;
	 if (dwperp=2) then begin
	    k:=k+1;
	    write(tx,ppar[k]:10:4,st[k]:2);
	 end;
	 writeln(tx,naoch[i,oche[i,j]]/natoms:10:4);
      end;
   end;
   writeln(tx,'end_dfit_params -------------------------------------------');
   writeln(tx,'start_ppos_params -----------------------------------------');
   export_params(tx);
   writeln(tx,'end_ppos_params -------------------------------------------');
end; { paradump }
procedure randomize_params(what	:integer);
Var
   p,i,j,o,l : integer;
   c	     : myreal;
begin
   {reset DW params to averages so they dont bias}
   if(what<>1) then begin
      for i:=1 to maxorb do begin
	 for j:=1 to noche[i] do begin
	    for l:=1 to nndw do begin
	       adwa[i,oche[i,j],l]:=adwa[i,oche[i,j],l]*(randu(seed)-0.5);
	       if(adwa[i,oche[i,j],l]<dwa0)then adwa[i,oche[i,j],l]:=2.0*randu(seed);
	       if(adwa[i,oche[i,j],l]>dwa1)then adwa[i,oche[i,j],l]:=dwa1-2.0*randu(seed);
	    end;
	 end;
      end;
   end;
   {random displacement in Angstrom by sigmc for each param}
   if (what>0) then begin
      for i:=1 to maxorb do begin
	 for o:=1 to 3 do if (o3pxyz[i,o]>0) then begin
	    p:=o3pxyz[i,o];
	    {normalize properly in-plane and z- components}
	    if(o<3)then c:=1/lpar[1] else c:=multi/lpar[2];
	    {loop through chemistrie}
	    for j:=1 to nn_pxyz[p] do
	       dpxyz[p,j]:=dpxyz[p,j]+sigmadr*c*(randu(seed)-0.5);
	 end;
      end;
   end;
end;
procedure dfit(xyzfit,dwfit : integer; Var epot,rf,rwf:myreal;Var chi3:v3);
var
   i,k,iter,m,n,p,o,nml,j,l : integer;
   chio,chi,ala,alao,cmp,c  : myreal;
   s1,s2,s3,sx		    : mystring;
   udw			    : boolean;
   se			    : rsites;
   st_save		    : gllista;
   mychi		    : myreal;
   move			    : integer;
begin
   OCH_STATUS;
   {xyzfit : 0=dont move; 1=move, no energy; 2=etarg mode 3=lambda mode}
   {dwfit: 0=keep fixed; 1=fit}
   move:=0; if(xyzfit>0) then begin
      move:=1;
      if(dwfit>0) then move:=2;
   end;
   {writeln('dfit: ',xyzfit:1,' ',dwfit:1,' ',move:1);}
   st_save:=st;
   udw:=((dwfit=2) and unidw);
   ORB2PARA(move);
   if(move=1) then begin
      m:=psca-1; if (dwperp=1) then m:=psca-2;
      for i:=pdw to m do st[i]:=0;
   end;
   if udw then unidw_setup(ndwdat) else ndwdat:=0;
   nedatap:=0;
   if(move>0) then begin
      if (xyzfit=2) then begin
	 nedatap:=1;
	 fo[ndat+ndwdat+1]:=etarg;
	 fos[ndat+ndwdat+1]:=lambdafit/sqrt(ndat*chinorm);
      end;
      if(xyzfit=3) then begin
	 efit_setup(ndwdat,nedatap);
      end;
   end;
   nml:=0;
   iter:=0; ala:=-0.1; alao:=ala; 
   repeat
      nml:=nml+1;
      npfix:=0;
      MAKE_PARA_LIST(st);
      repeat
	 iter:=iter+1;
	 if udw then PARA2AVEDW;
	 MRQMIN(q,fo,fos,ndat+ndwdat+nedatap,ppar,ma,plist,mfit,covar,alpha,ma,chi,ala);
	 if (abs(chio-chi)<dchi*ndat) then n:=n+1 else n:=0;
	 alao:=ala;
	 chio:=chi;
	 {writeln('mrq ',iter:1,' ',ala:1:5,' ',chi/ndat:1:2,' ',n:1);}
      until ((iter=lvmiters) or (n>10));
      {writeln('mrq inl(',nml:1,'): chi=',chi/ndat:1:3);}
      for i:=1 to mfit do if (st_fix[i]=1) then st[i]:=0;
   until npfix=0;
   {this sets displacements from ppar}
   if(move>0)then force_ft_init:=true;
   for i:=1 to ndat+ndwdat+nedatap do funcs_mrq(i,ppar,fc[i],ppard,0);
   chicomp(chi3);
   rwfac(fo,fc,fos,ndat,rf,rwf,mychi);
   if (xyzfit=2) then begin
      epot:=fc[ndat+ndwdat+1]*natoms;
   end else UPPAW(cueff,epot,se);
   PARA2ORB;
   st:=st_save;
end; { dfit }
procedure dfit_oneDW(orb,ch: integer; Var rf,rwf,mychi:myreal);
var
   i,k,l,iter,m,n,p,pch	: integer;
   chio,chi,ala,alao	: myreal;
   st_save		: gllista;
begin
   st_save:=st;
   for i:=pdw to psca-1 do st[i]:=0;
   pch:=poc(orb,ch);
   if unidw then begin
      k:=0; {dw-datapoint counter}
      for l:=1 to nndw do begin
	 if(peledw[ch]>0) then begin
	    p:=o2p[i,pch,l];
	    if(p>0) then begin
	       st[p]:=1;
	       k:=k+1;
	       pou[k]:=o2p[orb,pch,l];
	       fo[ndat+k]:=adwglo[ch,l];
	       fos[ndat+k]:=dwsig[ch];
	    end;
	 end;
      end;
      ndwdat:=k;
   end;
   iter:=0; ala:=-0.1; alao:=ala; 
   MAKE_PARA_LIST(st);
   repeat
      iter:=iter+1;
      MRQMIN(q,fo,fos,ndat+ndwdat,ppar,ma,plist,mfit,covar,alpha,ma,chi,ala);
      if (abs(chio-chi)<dchi*ndat) then n:=n+1 else n:=0;
      alao:=ala;
      chio:=chi;
      {writeln('mrq ',iter:1,' ',ala:1:5,' ',chi/ndat:1:2,' ',n:1);}
   until ((iter=lvmiters) or (n>10));
   for i:=1 to ndat+ndwdat do funcs_mrq(i,ppar,fc[i],ppard,0);
   rwfac(fo,fc,fos,ndat,rf,rwf,mychi);
   st:=st_save;
end; { dfit_oneDW }
procedure dwread(iname : integer);
var
   i,j,k,l,sta,i0,o0,po	  : integer;
   txt			  : text;
   rf,rwf,ac,pr		  : myreal;
   t_dwa,t_dwp,gdwp,ndwa0 : integer;
begin
   {Transfers dw parameters from one approximant to another,
   BUT: binding MUST be the same}
   if (iname=1) then reset(txt,'INPUT/DFIT/ini') else reset(txt,'INPUT/DFIT/ppar.dump');
   readln(txt,scal,st[psca],aco,t_dwa,t_dwp,dwpg);
   ndwa0:=1;
   if (t_dwa>0) then ndwa0:=2;
   if (t_dwa=3) then ndwa0:=3;
   ppar[psca]:=scal;
   if (aconst=1) then ppar[psca+1]:=ac;
   if (dwperp=1) then ppar[psca-1]:=dwpg;
   repeat
      read(txt,i0,o0);
      for l:=1 to ndwa0 do begin
	 read(txt,pr,sta);
	 adwa[i0,o0,l]:=pr;
	 iadwa[i0,o0,l]:=sta;
      end;
      if (t_dwp=2) then begin
	 read(txt,pr,sta);
	 adwp[i0,o0]:=pr;
	 iadwp[i0,o0]:=sta;
      end;
      readln(txt);
      write('DWR ',i0:4,o0:4,adwa[i0,o0,1]:1:4,' ',adwa[i0,o0,2]:1:4,' ',adwa[i0,o0,3]:1:4);
   until eof(txt);
   close(txt);
end; { dwread }
procedure dfit_glo;
Var
   i,j,k,m,ch : integer;
   fmax,epot  : myreal;
   chi,rf,rwf : myreal;
   se	      : rsites;
Var	      
   iter,n,p	 : integer;
   chio,ala,alao : myreal;
begin
   glodw:=true;
   {ignore load_dw}
   for i:=1 to atps do for j:=1 to 2 do adwglo[i,j]:=0.0;
   GLO2PARA;
   m:=0;
   fmax:=0.0;
   force_ft_init:=true;
   for i:=1 to ndat do begin
      FUNCS_MRQ_GLO(i,ppar,fc[i],ppard,m);
      if (fc[i]>fmax) then fmax:=fc[i];
      {writeln(i:1,' ',fc[i]:1:3,' ',fmax:1:3);}
   end;
   ppar[psca]:=1/fmax;
   iter:=0; ala:=-0.1; alao:=ala; 
   repeat
      iter:=iter+1;
      MRQMIN(q,fo,fos,ndat,ppar,ma,plist,mfit,covar,alpha,ma,chi,ala);
      if (abs(chio-chi)<dchi*ndat) then n:=n+1 else n:=0;
      alao:=ala;
      chio:=chi;
      {writeln('GLO: ',iter:2,' ',ala:1:5,' ',chi/ndat:1:2);}
   until ((iter=lvmiters) or (n>10));
   for i:=1 to ndat do funcs_mrq(i,ppar,fc[i],ppard,0);
   PARA2GLO;
   rwfac(fo,fc,fos,ndat,rf,rwf,chi);
   writeln('GLO_INIT: ',epot/natoms:1:5,' ',ndat:1,'  ',chi:1:3,' ',rwf:1:4,' ',rf:1:4);
   write('DW factors: ');
   for i:=1 to nnele do begin
      write(ele[i]:1,' ',adwglo[ele[i],1]:1:3,' ');
      if (aniso>0) then write(adwglo[ele[i],2]:1:3,'  ') else write('  ');
   end;
   writeln;
   glodw:=false;
end;
procedure dfit_setup;
Var
   i,j,k,m,ch	   : integer;
   fmax,epot,elas  : myreal;
   chi,rf,rwf,chix : myreal;
   rfint,rwfint	   : myreal;
   se,sel	   : rsites;
   chi3		   : v3;
begin
   OCH_STATUS; {initialize naoch and nnat and natoms}
   if(glodw) then DFIT_GLO;
   glodw:=false;
   if(load_dw>0) then DWREAD(load_dw);
   {initialize fourier transform, Debye Wallers...}
   ORB2PARA(0);
   ndwdat:=0;
   nedatap:=0;
   force_ft_init:=true;
   for i:=1 to ndat do FUNCS_MRQ(i,ppar,fc[i],ppard,0);
   {Note: this resets scale factor unless dwinp}
   if not(load_dw=1) then begin
      fmax:=0.0;
      for i:=1 to ndat do begin
	 if (fc[i]>fmax) then fmax:=fc[i];
      end;
      ppar[psca]:=1/fmax;
      ORB2PARA(0);
      DFIT(0,1,epot,rf,rwf,chi3);
   end;
   if(setapsig=1) then begin
      INTERNAL_RWFAC(rfint,rwfint,chix); {this sets fcs for approximant}
      for i:=1 to ndat do fos[i]:=fosd[i]+apsig*fcs[i];
      DFIT(0,1,epot,rf,rwf,chi3); {redo the fit}
   end;
   writeln('INIT: ',ndat:1,'  ',epot/natoms:1:6,' ',chi3[1]:1:3,' ',rwf:1:4,' ',rf:1:4,' # ndat E chi wR R');
   {internal_rwfac(rf,rwf,chi);
    writeln('                      ',chi:1:3,' ',rwf:1:4,' ',rf:1:4,' # internal chi wR R : ');
    elastic_E(cueff,epot,elas,se,sel);
    writeln('E_elas=',elas/natoms:1:6,' eV/at -> estim E_min=',(epot-elas)/natoms:1:6,' eV/at');}
   if(unidw) then unidw_setup(ndwdat);
   PARA2AVEDW; {regardless of the unidw, just get the averages}
   if (debug[13] and (unidw)) then begin
      write('dfit_setup: biasing DW uniformity : ',ndwdat:1,' extra datapoints (chem:');
      for i:=1 to neledw do write(' ',eledw[i]:1);
      writeln(')');
   end;
   PARA2ORB;
end;
procedure dfit_init;
var
   i,j,l : integer;
begin
   READ_CROMER;
   if (debug[10]) then write('Reading INPUT/DFIT/dat.hkl: ');
   READ_DATA(qplim,sigmalim,flim); {diffraction data}
   {chinorm:=0; for i:=1 to ndat do chinorm:=chinorm+1/fos[i]/fos[i];}
   chinorm:=1;
   {writeln('dfit_init: chinorm=',chinorm:1:5);}
   {init dw}
   for i:=1 to maxorb do
      for j:=1 to z_ele do begin
	 for l:=1 to nndw do begin
	    adwa[i,j,l]:=randu(seed);
	    iadwa[i,j,l]:=1;
	 end;
	 adwp[i,j]:=0;
         iadwp[i,j]:=1;
      end;
   DFIT_SETUP;
end; { dfit_init }
procedure hkl_out(fn :mystring);
Var
   i,o : integer;
   txt : text;
begin
   rewrite(txt,fn);
   writeln(txt,'# Fo Fc sig_used sig_obs sig_app Q Qp dQ po');
   for i:=1 to ndat do begin
      write(txt,fo[i]:1:5,' ',fc[i]:1:5,' ',fos[i]:1:5,' ',fosd[i]:1:5,' ',fcs[i]:1:5);
      write(txt,' ',i:1,'  ');
      for o:=1 to 3 do write(txt,hkl[i,o]:1,' ');
      writeln(txt,' ',q[i]:1:5,' ',qpe[i]:1:5,' ',dq[i]:1:6,' ',phkl[i]:1);
   end;
   close(txt);
end; { hkl_out }
procedure dcmc(nsteps,kk : integer; Var epot,ueff,rf,rwf,mychi:myreal);
{ASSUMES BNN table is initialized!}
{displace only one orbit at a time}
Var
   i,j,k,l,o,il,norb	: integer;
   ep1,eini		: myreal;
   pchem,porb,nrect,mpa	: integer;
   cf,cmp,sg2,cs2,cs	: myreal;
   du0,x,y		: myreal;
   v,w,w1,drt,dpr	: v3;
   br			: v33;
   se			: rsites;
   accept		: boolean;
   nacc,lab		: integer;
   m1,m2,nao,naao	: integer;
   fts			: array[1..ndatap] of v2;
   ffs			: array[1..ndatap,1..t_orb,1..t_ele] of myreal;
   chi1,chini,chi	: myreal;
   fcsave		: glndata;
   lmbd,sigma		: myreal;
   c1,c2,cn2,echi2	: myreal;
   s1,s2,s3,sx		: mystring;
   pparsave		: glmma;
   xwt,rate,rmsd	: myreal;
   chi3			: v3;
   procedure update_par(m1,m2,isg: integer;dr:v3);
   Var
      o,k,pp : integer;
   begin     
      for o:=1 to 3 do if (n_para[m1,m2,o]>0) then begin
	 pp:=n_para[m1,m2,o];
	 for k:=1 to nn_pxyz[pp] do if (n_a[m1,m2,1]=ch_pxyz[pp,k]) then begin
	    dpxyz[pp,k]:=dpxyz[pp,k]+isg*dr[o];
	 end;
      end;
   end; { update_par }
   procedure ff_init(a: glmma);
   Var
      o,i,j,po,k     : integer;
      dwex,dwexp,cro : array[1..t_orb,1..t_ele] of myreal;
      dwexpg	     : myreal;
   begin
   {intensity and partial derivatives}
      for po:=1 to ndat do begin
	 dwexpg:=1.0;
	 if (dwperp=1) then begin
	    dwexpg:=exp(a[psca-1]*qpe2[po]);
	 end;
	 k:=pdw-1;
	 for i:=1 to maxorb do begin
	    for j:=1 to noche[i] do begin
	       if (nndw=1) then begin
		  k:=k+1;
		  dwex[i,j]:=exp(a[k]*ex[po]);
	       end;
	       if (nndw=2) then begin
		  k:=k+2;
		  dwex[i,j]:=exp(a[k-1]*qd2[po]+a[k]*qp2[po]);
	       end;
	       if (nndw=3) then begin
		  k:=k+3;
		  dwex[i,j]:=exp(a[k-2]*qd2[po]+a[k-1]*qp2[po]+a[k]*qm2[po])
	       end;
	       dwexp[i,j]:=1.0;
	       if (dwperp=2) then begin
		  k:=k+1;
		  dwexp[i,j]:=exp(a[k]*qpe2[po]);
	       end; 
	       cro[i,j]:=cromer(oche[i,j],q[po]/2.0);
	       ff[po,i,j]:=dwexpg*dwex[i,j]*dwexp[i,j]*cro[i,j];
	    end;
	 end;
	 if(dwperp=1)then k:=k+1;
	 if(k+1<>psca)then begin
	    writeln('ff_init : bad parameter map, k+1 psca -> ',(k+1):1,' ',psca:1);
	    halt;
	 end;
      end;
   end; { ff_init }
   procedure ft_upd(orb,ch,dupmod :integer);
   Var
      i,j	   : integer;
      s		   : i3;
      na,pch,po	   : integer;
      co,si,arg,ac : myreal;
      y,pi2	   : myreal;
      x		   : v2;
   begin
      pi2:=8.0*arctan(1.0);
      pch:=poc(orb,ch);
      if (aconst=1) then ac:=ppar[psca+1] else ac:=0.0;
      for po:=1 to ndat do begin
	 fcsave[po]:=fc[po];
	 for o:=1 to 2 do fts[po,o]:=ft[po,orb,pch,o]; {save}
	 for o:=1 to 2 do ft[po,orb,pch,o]:=0.0;
	 na:=0;
	 for i:=1 to tnodes do begin
	    for j:=1 to nan[i] do begin
	       if(n_a[i,j,1]>0) then na:=na+1;
	       if ((n_a[i,j,2]=orb) and (n_a[i,j,1]=ch)) then begin
		  arg:=0.0;
		  for o:=1 to 3 do arg:=arg+hkl[po,o]*(rt[i,o]+n_r[i,j,o]);
		  co:=cos(pi2*arg);
		  si:=sin(pi2*arg);
		  ft[po,orb,pch,1]:=ft[po,orb,pch,1]+co;
		  ft[po,orb,pch,2]:=ft[po,orb,pch,2]+si;
	       end;
	    end;
	 end;
	 for o:=1 to 2 do ft[po,orb,pch,o]:=ft[po,orb,pch,o]/na;
	 if(dupmod=0) then begin
	    for o:=1 to 2 do x[o]:=0;
	    for i:=1 to maxorb do begin
	       for j:=1 to noche[i] do if (aino(i,j)>0) then begin
		  for o:=1 to 2 do x[o]:=x[o]+ft[po,i,j,o]*ff[po,i,j];
	       end;
	    end;
	 end;
	 y:=sqrt(abs(x[1]*x[1]+x[2]*x[2]+ac));
	 fc[po]:=ppar[psca]*y;
      end; {end ddataloop}
   end; { ft_upd }
begin
   xwt:=cueff;
   lmbd:=lambdafit/sqrt(ndat*chinorm);
   sigma:=sigmadr;
   lab:=10;
   mpa:=psca; if (aconst>0) then mpa:=mpa+1;
   DFIT(0,2,epot,rf,rwf,chi3);
   chini:=chi3[1];
   UPPAV(eini,se); {update energy after the last step}
   ff_init(ppar);
   nacc:=0;rate:=0.0;naao:=0;
   sg2:=2*sigma*sigma;
   cs:=0.0; cs2:=0.0;
   echi2:=0.0;
   UPPAW(xwt,epot,se);
   writeln('#STARTLOOP ',eini/natoms:1:6,' ',epot/natoms:1:6,' ',chini:1:3,' ',rwf:1:4,' ',rf:1:4,' ',lmbd:1:3,' ',xwt:1:2,' ',sigma:1:2,' ',kk:1,' ',beta[kk]:1:1,' # E UEFF CHI RWF RF LAMBDA CUEFF LOOP BETA');
   {writeln(txout,'#STARTLOOP ',eini/natoms:1:6,' ',epot/natoms:1:6,' ',chini:1:3,' ',rwf:1:4,' ',rf:1:4,' ',natoms:1,' ',ndat:1,' ',lmbd:1:6,' ',xwt:1:3,' ',sigma:1:3,' ',kk:1,' ',beta[kk]:1:3,' # E UEFF CHI RWF RF NAT NDAT LAMBDA CUEFF LOOP BETA');}
   chi1:=chini;
   {writeln('CHINORM ',chinorm:1:4);}
   nrect:=round(rectime*nsteps*natoms);
   for il:=1 to nsteps*natoms do begin
      repeat
	 cmp:=randu(seed)*tnodes; {pick randomly tiling node}
	 m1:=round(cmp+0.5); 
	 cmp:=randu(seed)*nan[m1]; {pick randomly an atom}
	 m2:=round(cmp+0.5);
      until (n_a[m1,m2,1]>0);
      pchem:=n_a[m1,m2,1];
      porb:=n_a[m1,m2,2];
      nao:=naoch[porb,pchem];
      for o:=1 to 3 do drt[o]:=randu(seed)-0.5;
      repeat {gaussian random length (0,drmax)}
	 cf:=randu(seed)*sigma*5; {should be big enough}
	 y:=randu(seed);
      until (exp(-cf*cf/sigma/sigma/2)>y);
      cmp:=absv3(drt);
      for o:=1 to 3 do drt[o]:=drt[o]*cf/cmp;
      get_atom_basis(n_para[m1,m2],br);
      mxle3(br,drt,dpr); {parameters basis}
      {for o:=1 to 3 do dpr[o]:=0.0;}
      update_par(m1,m2,1,dpr); {forward}
      {energy}
      UPPAW(xwt,ep1,se); {loads new n_r and calculates energy}
      FT_UPD(porb,pchem,dwupdmode); {saves fc and ft}
      if (dwupdmode=0) then begin
	 rwfac(fo,fc,fos,ndat,rf,rwf,chi1);
      end;
      if (dwupdmode=1) then begin
	 for i:=pdw to mpa do pparsave[i]:=ppar[i]; {save DW+ params}
	 DFIT_ONEDW(porb,pchem,rf,rwf,chi1);
      end;
      if (dwupdmode=2) then begin
	 for i:=pdw to mpa do pparsave[i]:=ppar[i]; {save DW+ params}
	 DFIT(0,2,epot,rf,rwf,chi3);
	 chi1:=chi3[1];
      end;
      {xw:=(1+(natoms/maxorb-1)*xwt);}
      c1:=ep1-epot;
      c2:=lmbd*(chi1-chi);
      cn2:=c1*c1+c2*c2;
      du0:=(c1+c2)*beta[kk];
      accept:=true;
      if(du0>0)then begin
	 x:=exp(-du0);
	 y:=randu(seed);
	 if (y>x) then begin {reject}
	    update_par(m1,m2,-1,dpr); {backward}
	    accept:=false;
	    for i:=1 to ndat do begin
	       if (dwupdmode=0) then fc[i]:=fcsave[i];
	       k:=poc(porb,pchem);
	       ft[i,porb,k]:=fts[i];
	    end;
	    if (dwupdmode>0) then for i:=pdw to mpa do ppar[i]:=pparsave[i];
	 end;
      end;
      if(accept)then begin
	 {if (debug[8]) then
	  writeln(il:2,'   ',ep1/natoms:1:6,' ',(ep1-epot):1:4);
	  if ((il mod save_e)=0) then begin
	  write(recu,(ep1-epot):1:6,' ');
	  writeln(recu,'   ',kk:1,' ',lab:1);
	 end;}
	 epot:=ep1;
	 chi:=chi1;
	 if(il>nrect) then begin
	    nacc:=nacc+1;
	    cs2:=cs2+cf*cf*nao;
	    cmp:=(c1*c1)/cn2;
	    if(cn2>1e-10)then begin
	       echi2:=echi2+cmp;
	    end else echi2:=echi2+0.5;
	    naao:=naao+nao;
	    {writeln(il:1,' ',c1*beta[kk]:1:4,' ',c2*beta[kk]:1:4,'   ',cmp:1:4,' ',rwf*100:1:4,' ',beta[kk]:1:2);}
	    {writeln(txout,(epot):1:4,' ',rwf*100:1:4,' ',kk:1,' ',cmp:1:4);}
	 end;
	 {mystat(kk,lab,epot);
	  if(epot<emin)then begin
	  export_config_if(emin,epot/natoms,beta[kk],sname[9]);
	 end;}
      end;
   end; {main nsteps loop}
   ueff:=epot;
   rmsd:=sqrt(cs2/naao);
   UPPAV(epot,se); {update energy after the last step}
   {if debug[8] then begin}
   if(naao>0)then writeln('rms-displacement: ',sqrt(cs2/naao):1:4);
   {end;}
   rwfac(fo,fc,fos,ndat,rf,rwf,chi);
   xyz2ppar(ppar);
   {force_ft_init:=true;
    for j:=1 to ndat do funcs_mrq(j,ppar,fc[j],ppard,0);
   rwfac(fo,fc,fos,ndat,rf,rwf,chi);}
   rate:=nacc/(nsteps*natoms*rectime);
   writeln('#ENDLOOP: ',epot/natoms:1:6,' ',ueff/natoms:1:6,' ',rate:1:4,' ',sqrt(cs2/naao):1:4,' ',echi2/nacc:1:4,' ',chi:1:2,' ',rwf:1:4,' ',rf:1:4,' # E UEFF RATE RMSD EWT CHI RWF RF');
   {writeln(txout,'#ENDLOOP: ',epot/natoms:1:6,' ',ueff:1:6,' ',rate:1:4,' ',sqrt(cs2/naao):1:4,' ',echi2/nacc:1:4,' ',chi:1:2,' ',rwf:1:4,' ',rf:1:4,' # E UEFF RATE RMSD EWT CHI RWF RF');}
   mychi:=chi;
end;

end.

